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Abstract 

In the focusing problem we study a solution of the porous medium 
equation ut = A (n™') whose initial distribution is positive in the ex- 
terior of a closed non-circular two dimensional region, and zero inside. 
We implement a numerical scheme that renormalizes the solution each 
time that the average size of the empty region reduces by a half. The 
initial condition is a function with circular level sets distorted with a 
small sinusoidal perturbation of wave number k > 3. We find that 
for nonlinearity exponents m smaller than a critical value which de- 
pends on k, the solution tends to a self-similar regime, characterized by 
rounded polygonal interfaces and similarity exponents that depend on 
m and on the discrete rotational symmetry number k. For m greater 
than the critical value, the final form of the interface is circular. 
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1 Introduction 



This work deals with solutions to the nonlinear diffusion equation 

l-AK"), (1) 

where m > 1 is constant and A denotes the Laplace operator in R"'. Equa- 
tion (|I]) arises in various problems, such as the spreading of viscous gravity 
currents ||, the diffusion of strong thermal waves and the isentropic 
flow of an ideal gas in a homogeneous porous medium In the latter 

application, u represents the scaled density of the gas. For m = 1 equation 
([^) is the classical heat conduction equation and there are other applications 
which involve values of m < 1. One of the salient features of the "slow diffu- 
sion" case m > 1 is the occurrence of interfaces moving with finite velocity, 
that separate empty regions where m = from regions where m > 0. More 
detailed information about the properties of solutions to equation ([l|) can be 
found in 

In order to discuss the behavior of solutions to equation (|l|) it is convenient 
to replace the independent variable u with 

^ m-l 

V = u . 



m — L 

In the gas flow interpretation, v represents the scaled pressure via the ideal 
gas law. If u satisfies (P then v satisfies 

= (m-l)vAv+\Vvf , (2) 

ot 

where V is the gradient operator in R'^. 

Here we are interested in the focusing problem in which we solve equation 
(0) starting from an initial distribution which is positive in the exterior of 
a closed bounded region D and zero inside D (see Fig. p. At some finite 
time T > 0, called the focusing time, v will for the first time be positive 
throughout the initially empty region D. 

In IP it is shown that for each m G (l,C)o) there is a unique similarity 
exponent (3{m,d) G Q, ij and a one-parameter family of functions Fc such 
that the functions 

v^{x,t) = {T-tY^-'Fjj^}^y (3) 
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with 13 = (5 {m,d), are self-similar solutions to (Q). Moreover, there exists a 
7((i, m) such that 



= for < e < (^) 
> for e > (^Y 



/3 



where 7 = 'j{d,m). Thus the Wc(a;, t) are focusing solutions to equation 
with interfaces given by 

\x\ 




(T - 1)^ 

The focusing problem is well studied in the radially symmetric case. In 
particular, it is proven in Q that some member of the Graveleau family 
(|^) describes locally, to leading order, the focusing behavior of essentially any 
radially symmetric focusing solution to the pressure equation (Q). The situ- 
ation is much more complicated in the absence of radial symmetry. Physical 
experiments involving convergent gravity flows (m = 4, c? = 2) followed by 
numerical experiments indicate that both large and small deviations from 
rotational symmetry may be amplified as the flow tries to fill the hole. A 
formal linear stability analysis for general m and d (cf. and Appendix 2) 
shows that the Graveleau solutions are indeed unstable, at least when m is 
close to unity, and that the number of unstable modes increases as m \ 1. 
This suggests that a sequence of bifurcations occurs as m decreases from 00 
to 1. The existence of these bifurcations is proved in [Q]. The bifurcating 
solutions are self-similar, but not axially symmetric. They are invariant with 
respect to the action of the group 0{d — l, R) of {d—l)x (d — l) real orthogo- 
nal matrices. In this paper we restrict our attention to the plane case d = 2. 
In this case the results of show that for any rh G (1, 00) there exists an 
integer k^{m) G (2, 00) such that a symmetry breaking bifurcation from the 
radially symmetric Graveleau solutions occurs at some G {l^fh) for each 
k > k^:{fn). The bifurcating solutions are non-radial self-similar focusing so- 
lutions possessing A;-fold symmetry, i.e., the symmetry of cos{k6). Numerical 
studies show that for each k > 3 the bifurcation point m = irtk is unique so 
that the Graveleau solutions are linearly stable with respect to perturbations 
with A;-fold symmetry for m > and unstable for > m > 1. Moreover, 
the ruk are ordered 

00 > rrij, > > ... > > ... \ 1. 
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Crude estimates for the first four rrik are given in Table I (See also Table IV 
in Appendix 2). 

k estimate 

3 mg e (1.69, 1.7) 

4 m4 G (1.32, 1.321) 

5 ms G (1.18, 1.19) 

6 me G (1.12, 1.13) 

Table I: Estimates of the values of for which the circular self-similar 
solution bifurcates into a non-circular solution. 

The results in ^ give no information about stability with respect to 
perturbations with wave number k = 2. The experimental and numerical re- 
sults in suggest instability. This is confirmed by numerical linear stability 
analysis which shows that the Graveleau interface is unstable with respect 
to perturbations with wave number 2 for all values of m. In particular, there 
appear to be no bifurcating branches of self-similar focusing solution with 
the 2-fold symmetry. We discuss this briefly in Section 5 of this paper and 



in more detail in reference ITB 



In this paper we carry out detailed numerical studies of focusing solutions 
to equation (H) whose initial distributions have interfaces which are circles 
with small perturbations. Mainly, we deal with initial conditions whose in- 
terfaces are of the form 

r = a(l + £:cos(A;6')) (4) 

with e -C a, although in Section 5 we will consider perturbations with mixed 
modes. By using a numerical renormalization technique inspired by the pi- 
oneering work of Chen and Goldenfeld |TT| we are able to follow the 
evolution of the interface to times very close to the focusing time. Thus we 
are able to obtain very detailed information about the asymptotic form of the 
solution as it focuses. As we shall see below, for single mode perturbations 
with fc-fold sjTiimetry, the numerical results indicate that the leading term 
in the focusing asymptotics is a self-similar solution of the form 

v = iT-tr-W.^-^^,9^ (5) 

where c is a parameter, 6 = 6 {k^ m) is the similarity exponent, and Vc satisfies 

K(C, 0) = K(C, ^ + -^) forn = 1, A: - 1. (6) 
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Moreover, the focusing interface is asymptotically of the form 

r = (T - tfA {6) , 

and S{k,mic) = /3{2,mk), where f3 is the similarity exponent for the Grav- 
eleau solution and nik is the bifurcation point found in 0. Therefore we 
conclude that the functions given by (|^) are the bifurcating solutions found 
in i. 



2 Numerical scheme and renormalization pro- 
cedure 

In order to solve Eq. (H) numerically, we discretize a circular domain D = 
[0,R] X [0,27r/A;] with a uniform polar grid of interval sizes 6r = R/Nr and 
66 = 271 /Nq. The numerical solution is stored in a matrix as v{r,6,t) = Vij 
with r = i6r and 6 = j66. In order to integrate the PDE in time, we use an 
explicit Euler scheme 

v,j{t + St) = Vij{t) + StH.j (7) 

where Hij is a finite differences approximation of the derivatives of the right 
hand side of Eq. (||). We compute the derivatives of the term |Vf |^ = f ^ + 
vj/r"^ with a second order upwind ENO (essentialy non-oscillatory) scheme 



1^ , p!5| (see Appendix 1). This method guarantees that the numerical scheme 
will be able to describe accurately the discontinuities on the first derivative 
that spontaneously appear in the case m = 1 (Hamilton- Jacobi limit) and at 
the interface. Due to the diffusive nature of the Laplacian term in Eq. (^, 
the corresponding derivatives are computed with standard centered second 
order approximations, 

Av = Vrr-\ \- ^ - (8) 

(5r2 ^ 26rri 69^rf ' ^ ' 

Near the interface there is a discontinuity in the first derivative that may 
generate numerical errors when we compute the Laplacian with Eq. (^. 
Therefore, for the grid points Vij situated at a distance smaller than 25r 
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from the interface, the Laplacian is computed by hnearly extrapolating in 
the variable r the values of the Laplacian from the nodes of behind, which 
were previously computed with Eq. (^. 

At the boundaries 6 = ±7r /k we apply the boundary condition of period- 
icity v{—Tr/k) = v{7T/k). This is simply a convenience, since computations 
without this forced symmetry yield equivalent results. At the boundary 
r = Rwe apply the boundary condition Vrr = 0. This is equivalent to a first 
order linear extrapolation of the ghost points vn^+ij outside R, which are 
needed in order to compute the derivatives at R. Finally, at r = we set 
v = 0. 

We start the integration with a rather arbitrary initial condition. 



V 



^^\r,e,to) =r -ao{l + ecoske) , e = 0.1 (10) 



which describes a function whose contour lines are perturbed circles, and 
which interface is given by Eq. @ . The exact form of the initial condition is 
not very critical, because the asymptotic solution only depends on k (which 
determines the symmetry) and m, as was verified numerically. 

We integrate the diffusion equation over a sequence of time intervals 
{tn,tn+i) starting with n = 0, and renormalize the solution at each right 
hand end point t = before continuing the integration. The renormal- 
ization times t^+i are taken to be the times when the average radius of the 
interface]^ ^ 

a(t) = 7^ / " a{e,t)de (11) 



27r Jo 

reaches half of its initial value. The renormalized solution is defined to be 

v^^^'\r,e,t^+,) = ■t;(")(r/2,^,t„+i) (12) 

for n = 0, 1, 2, .... This transformation is performed by linearly interpolating 
the values of the grid, 

vlj^^^ = Z^"^^ ■ vljlj for even i 



and 



vt'^ = Z^-"^ ■ + v^l,^^)/2 for odd ^. 



^ The interface a{9,t) is defined as v{a{d,t),6,t) = 0. fn order to avoid numerical 
problems due to the numerical diffusion near the interface, we extrapolate the positive 
values of w up to w = 0. 
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The constant Z^"^^ is taken as the reciprocal of the maximum value of the 
function v in the renormalized domain of integration 

Z(")= fmaxt;(")(r,^,t„+i)^ . (13) 

\r<R/2 J 

The superscript (n) indicates how many renormalizations we have made up 
to the time tn+i- The errors introduced in the linear interpolation are of 
second order, the same as in the discretization of the derivatives. It is very 
important to determine the time tn+i accurately, and this is done by first 
detecting the exact time when the interface crosses half of its initial value at 
the previous renormalization, and then using this information to re-compute 
the time step 6t and the solution such that the interface reaches exactly the 
half of the initial value. 

We summarize the procedure as follows: 

1. Initialize f (°)(r, 6*, to) with Eq. (|10D. Let n = 0. 

2. Starting from t = tn solve Eq. (^) with standard finite differences until 
the time t = tn+i where a(t„+i) = a{tn)/2. Here we have v^"\r, 9, t„+i). 

3. Renormalize w("+i)(r, ^, t„+i) = ■ t;(")(r/2, ^, t„+i). 

4. While 

(14) 

let n = n + 1 and return to step 2. 

where r is a tolerance. In the experiments reported below, we have set 
r = 10"^ 

We found that when we repeat this procedure n times, the solution typi- 
cally converges for n > 50. This convergence is a necesary condition for the 
existence of a self-similar solution. In Fig. (^ we show the convergence of 
the iterative procedure by plotting the difference of succesive approximations 
given by Eq. (|14[) and the departure of Z*^") from its asymptotic value. 

In order to compute the exponent 6 we note that once the scheme has 
converged. 
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and 

(here we are ignoring the numerical errors). Then, from Eq. (|12D, 

v^''\r,9,Q = Z-v^''\r/2,e,tr,+,). (15) 

Thus, the solution at time tn+i is proportional to the solution at time t„ with 
the distances scaled by a factor two. 
By substituting v from Eq. (^, 

(here we have put T = 0) it is clear that 

tr. = 2'/\+i and {-t^f-' = Z{-t^+,f-\ (17) 

Solving these algebraic equations we find that the similarity exponent is given 
by 

6-' = log, (I) . (18) 

3 Validation 

1. In order to show the accuracy of this numerical computation, we com- 
pare in Table II the numerical solution for the circularly symmetric 
case with the Graveleau solution, and we find a good agreement in the 



8 



exponents {S = P in this case, see Eqs. ^ and |^) 0. 



1 
1 
1 
1 
1 
1 
1 
1 
1 
1 
2 



m 




1 

,2 
3 
,4 
,5 
6 
7 
,8 
9 




Numeric 
1.00000 
0.976749 
0.956533 
0.938768 
0.923066 
0.909103 
0.896595 
0.885311 
0.875052 
0.865675 
0.857052 



'Exact' 
1 



0.976744 
0.956517 
0.938756 
0.923071 
0.908983 
0.896346 
0.884940 
0.874546 
0.865052 
0.856333 



Table II: Exponents for the circular case, R = 1, Nj. = 300 

This comparison checks the accuracy of the renormalization procedure 
and the temporal integration. 

2. We also computed the solutions corresponding to the Hamilton- Jacobi 
case 171 = 1, and we compared the numerical result with the exact 
solution. The interfaces of the exact solution are regular polygons of k 
sides (A; > 3), and the function w is a set of k planes of the form 

f (r, e, t) = c{x + a{t)) = c(r cos 9 + a{t)) -Ti/k<e < njk (19) 

vir, 9 + 2'Kn/k, t) = v{r, 6, t) for n = 1 ..k — 1 

where a{t) = ct and c is a constant . The level lines are regular polygons 
and the similarity exponent is 5 = 1. This solution exhibits a set of 
k discontinuities in the derivatives at the apeses of the polygonal level 
sets. The level sets of the numerical solution are shown in Fig. ^ for 
A; = 3, 4, and 5. The numerical solution gives the value of the exponent 
with an absolute error of the order of 10~^ when a discretization of 

^We implemented the numerical scheme in FORTRAN language and we made the 
computations on an IBM RS/6000 computer 



9 



Nj. = 300 and Ng = 300 is used. The relative departure of the numerical 
solution from the theoretical solution is 



A, 



1 



1 



// 



t;(")(r,0,t„) -v{r,9,Qfr dr d9 



■n 



7r(i?/2)2 



v.. 



max 



r<R/2 



(20) 



where v is computed with Eq. ( |T9| ) and Vmax is the maximum value of 
the solution in the integration domain. The value of c is obtained from 
the numerical solution, c = f^"'^(a(0, t), 0, t). This difference is shown 
in Fig. (Q) as a function of the number of renormalizations n. Clearly, 
the numerical computation converges to the exact solution with relative 
errors smaller than 10~^. 

This comparison checks the accuracy of the two-dimensional evaluation 
of the {Vvf term and the renormalization routine. 

3. We also validated the code for the case of two colliding plane fronts 
{d = 1), where the exact solution is 



The exact value of the exponent is 1 for all m. For m = 2 the numerical 
exponent are close to 1.0024, and the level lines depart from the straight 
theoretical lines in less than one part in 500. In this case, the Laplacian 
term is zero, however, since the numerical code is written in polar 
coordinates, it serves as a test for the evaluation of the terms in the 
Laplacian. 

4 Non-circular self-similar collapse 

We summarize here the results for non-circular collapses, starting with the 
initial condition Eq. ([lO|) . The numerical solution shows that the final shape 
of the interfaces are rounded polygons. The shape of the interfaces may be 
characterized by the ratio between the maximum to minimum radii of the 
interface 



v{x, t) = X + t. 




For instance, for a polygon of k sides, I is equal to cosvr/fc. 
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The computations were performed in a domain of radius R = 1 with 
Oq = 0.1665 and e = 0.1. The number of grid points in the r direction 
was Nr = 300 and the effective number of points in the 6 interval (0,27r), 
Ng = 300. Table III shows the exponents that result from this computation. 
The value of / is indicated in parentheses. 



m k = 3 k = A k = 5 

1.0 1.0000 (0.5029) 1.0000 (0.709) 1.0000 (0.810) 

1.1 1.0000 (0.5029) 1.0000 (0.709) 0.9997 (0.811) 

1.2 1.0000 (0.5027) 0.9985 (0.710) 

1.3 1.0000 (0.5025) 0.9587 (0.770) 

1.4 0.9996 (0.5021) 

1.5 0.9964 (0.5027) 

1.6 0.9839 (0.5075) 

Table III: Exponents for the non-circular collapse 

The dashes indicate that for this values of k and m the solution tends 
to the circular case (Graveleau solution). This result serves as another vali- 
dation, since the solution becomes unstable for the values of m reported in 
Appendix 2 (see also Table I). Note that for m = 1.3 and k = 4, the aspect 
ratio I becomes closer to the circular case (= 1) because we are near the 
bifurcation. In Fig. (§) we plot the exponents for the radially symmetric 
case, and for the non-circular cases = 3, 4, 5 and 6. Note the continuity 
of the curves at the stability limits. The case k = 2 is not included because 
the corresponding evolution is not self-similar. In Fig. ^ we show typical 
asymptotic solutions for three different values of k. 

We observe that the form of the asymptotic solutions only depend on 
the value of k that determine the symmetry and m, even when the initial 
condition has a hole with a more complex shape than the described by Eq. 



5 Collapse without self-similarity 

The computations shows that self-similar solutions are obtained when the 
initial condition has rotational k-fold symmetry, with k > 2. When the 
initial condition does not satisfy this symmetry condition, then the collapse 
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may not be self-similar. In []^ the experiments indicate that when the initial 
shape is elongated in one direction, then there is not a self-similar solution. 
Moreover, the dimensions of the mayor axis and the minor axis seem to follow 
different power laws. 

In Fig. 0-a we consider the following example with initial condition 

y(o)(r,^,to) =r-ao(l -£cos3^-ecos5^), £ = 0.05 (21) 

Clearly, this shape does not have k-fold symmetry for any k > 2. This combi- 
nation of modes yields a slightly elongated figure. The evolution for m = 1.5 
leads to the intermediate shapes shown in Fig. ^(b)-(c). Initially the 3-fold 
component of the initial condition dominates and the shape is almost trian- 
gular, but, since the initial shape was slighty ellongated in the x direction, 
the final shape becomes increasingly ellongated. As far as the computation 
shows, the solution does not reach a self-similar regime in this case. The 
shape of the hole becomes increasingly elongated, until the numerical grid 
cannot accurately resolve the shape of the interface. We can interpret this 
result in terms of the interaction of the modes: the new mode that accounts 
for the ellongation, namely = 2, is created from the non-linear interaction 
of the modes 3 and 5. 

If instead of the modes 3 and 5 we start with the modes 3 and 6, the shape 
will have 3-fold symmetry, and a self-similar solution for the most dominant 
mode {k = 3) is eventually reached. Further discussions of these cases will 
appear in ref. |T^. 
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6 Appendix 1 

In order to compute the |Vf| term, we use the ENO scheme, which has 
been used very succesfuUy in the numerical solution of Hamilton- Jacobi 
equations |T^, |T^. The ENO scheme is an adaptive stencil interpolation pro- 
cedure which automatically obtains information from the locally smoothest 
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region, and hence yields a uniform high essentially nonoscillatory approxi- 
mation for picccwisc smooth functions. 

When computing the partial derivative respect to r corresponding to Vji 
at discrete nodes rj = hj, j = 0, ±1,±2, ... we first write the undivided 
differences 

'^{j-i — '^{j + 1, A; — 1) — w{j, k — 1) k — 1, d+ 1 

where d is the order of the approximation {d — 2 in our case). The ENO 
stencil-choosing procedure is optimally implemented by starting with — j 
and performing 

if \wiiij),k)\ > \w{i{j) - l,k)\ then = - 1 

for k = 2,...,d. Finally we compute the forward and the backward deriva- 
tives, 

1 

(dv/dr)^ ^ tYI - j): k)w{i{j), k) 

k=i 

and 

1 

(dv/dr) = tY. - 1) - j), k)w{i{j - 1), fc) 
fe=i 

where 

1 m+k—l m+k—1 

cM = - E n H). 

The same procedure is used to compute the derivatives in the 9 direction. In 
order to compute the upwind approximation for we write 



min(y{dv/dr) , -\- max {j^dv / dr)'^ 



Similar expressions are used to compute v^. This scheme is very useful when 
there are discontinuities in the first derivatives, as in our case at the interface 
and for m — 1. 



7 Appendix 2 

In this Appendix we describe the linearized stability analysis of the Gravclcau 
interfaces in dimension d — 2. Although the basic equations are derived and 



13 



analyzed in 0, there remain several questions which, at least for the present, 
can only be answered by numerical studies. We describe these studies here, 
and for the convenience of the reader give a brief outline of the derivation of 
the relevant equations. 

We begin by rewriting the porous medium pressure Eq. (|^) in self-similar 
coordinates. Let f (x, t) be a solution of Eq. and set 

^;(x,t) = (T-tf^-V(y,r), 

where 

and r = — ln(T — t). 



Then V satisfies the equation 
dV 

— = nVAV + \VVf -PyVV + {2(3 - l)V, (22) 

where n = m — 1. Let p = |y| and write the Graveleau solution (Eq. (|^)) 
normalized with c = 1 in the form 

t;i(x,t) = (r-tf^-V(p). 
Note that is a. steady state solution to Eq. (|22|) , i.e., 

i^>' + ^^V^'^ + {i^'f - (3pi^' + {213 - l)^ = 0. 

Moreover, 
and 



= for < p < 7"'' 
> for p > 7"^ 



iIj {p) ^ p as p — >• oo 

(cf. 0). We now restrict our attention to flows in two dimensions and 
introduce polar coordinates 

y = p (cos 9, sin 9) , 

where 

Vi + yi s-iid 9 = arctan — . 

yi 
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If we set W{p, 9, r) = V(j, r) the W satisfies 

dW 

- Pp-^ + (2/? - 1) W. (23) 
For p > the level curves of W are given implicitly by 

W{p,9,T)=p. 
Assuming that ^ 7^ we can solve for p to get 

p = S{p,e,T). 

All of the derivatives of W which appear in equation ( ]23|) can be calculated 
from the relation 

w{S{p,e,T),e,T)-p = o 

(cf. for details). Thus we derive the evolution equation for the level curves 
p = S: 

S^] -w- = np{S'^-S —] +— - — -2- 



dp J Or I dp"^ \dp J dp"^ \ / '^P dpd6 

(24) 

The Graveleau function ipip) is an increasing function for p > 7 ^ with 
range [0, 00). Thus we can invert to obtain p = \l/(p) which is an increasing 
function for p > with range [7"^, 00). Note that \1/ is a steady state solution 
to Eq. (p^) and satisfies 

np^^" = npi^'f - /3(^^')^ + + (2/3 - l)p^(^0^ (25) 

where ' = d/dp. Set S{p,9,t) = ^(p) + ^{p,9,t), where we assume that 
lel « 1- 
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To leading order ^ satisfies the linear equation 

+ {2^^3^' - ^2 _ 2^p^^' _ 3(2/3 - p- (26) 

Since the coefficients in Eq. ( p6D depend only on p there are solutions of the 
form 

e, t) = A{p) exp {tke + Ar) , 
where A satisfies the ordinary differential equation 

np^^A" + {2/3^'^3 - ^2 _ 2np^^' - 3(2/3 - ^' (27) 

+ {np{l - e) + {(3- A)*'} {^I'fA = 0. 

Here A; is the given wave number and A is an eigenvalue which must be 
determined. Since ^1/(0) = 7"^^ it follows from equation p4| that \l/'(0) = j (3- 
Thus, to leading order for p <^ 1, equation 12^ becomes 

npA' + A' + (/5 - \)-i'^<^A = 0. 

This equation has a regular singular point at p = 0, and possesses a unique 
analytic solution determined by the initial conditions 

A{0) = 1 and A'{0) = ^^7'^. (28) 

For p 1, we have 

^ 213 -Y 

Thus, to leading order for p ^ 1, equation |2^ becomes 

^ 2/5-1 (2/3-1)2 

Most solutions of this equation grow exponentially at infinity with 

In A ~ —p 2/3-1, 
16 



but there are also solutions with algebraic growth 



A = (^p^^ . (29) 

The eigenvalues are those values of A for which the solution of Eq. ( pTj) with 
initial values of Eq. (^8|) has the algebraic growth given by Eq. (|29|) at 
infinity. They are obtained numerically by a shooting technique. 

The eigenvalues of equation (^) are analyzed in detail in reference 0, 
where it is shown that for each m G (l,oo) they form a doubly infinite 
sequence Xkj{m). Here > is the wave number and j > is the number of 
zeros of the corresponding eigenfunction. Moreover, Aoi(m) = 0, Aoj(m) < 
for j > 2, Xkj{m) < for k > 1 and j > 1, and Afco(^) > for A; = or 
1. Bifurcations occur for those values of m where Xkj{m) changes sign, and 
this can only happen for j = and k > 2. In |^ it is shown that there are 
infinitely many bifurcations. Specifically, for every m' G (1, oo) there exists 
an integer k^,{m') > 2 such that for each k > k^{m') there is a bifurcation 
point rrik G (l,m'). The bifurcating solutions are not radially symmetric, 
but do possess the /c-fold symmetry of cos(fc6'). 

Numerical studies of the eigenvalues indicate that there are no bifurca- 
tions for wave number k = 2 for any value of m, and that there is a unique 
bifurcation point ruk for each k > 3. In particular. 



Afco("^) 
The ruk are ordered with 



> for m G (1, nik) 
< for m G (m^, oo) 



oo > ms > 777-4 >■■• > "mk > ■ ■ ■ \ 1. 



Table IV summarizes some of the numerical results for the sign of Xko{m). 
A minus sign in the position [k, m) indicates stability with respect to per- 
turbations with wave number k for the given value of m, while a plus sign 
indicates instability. The bifurcation points occur between adjacent m- 
values where the XkQ{rn) have opposite signs. Thus, for example, 7774 G 
(1.320,1.321). 
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m — 


1.12 


1.13 


1.18 


1.19 


1.320 


1.321 


1.69 


1.7 


2 


2 


+ 


+ 


+ 


+ 


+ 


+ 


+ 


+ 


+ 


3 


+ 


+ 


+ 


+ 


+ 


+ 


+ 






4 


+ 


+ 


+ 


+ 


+ 










5 


+ 


+ 


+ 














6 


+ 



















Table IV: Sign of Xko{m). Negative values indicate that the circular 
self-similar solution is stable for the corresponding pair of m and k. 
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v>0 




Figure 1: Sketch of a noncircular convergent flow. The fluid hes outside the 
interface and flows inwards. 
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Figure 2: Convergence of the numerical solution: the difference between 
succesive approximations tend to zero and the renormalization constant Z 
tends to a constant as the number of renormalizations n increases. In this 
case m = 1.5, Nr = 300, Ng = 300 and the initial shape is a perturbated 
circle with k = 3. 




Figure 3: Polygonal contour lines of the asymptotic solutions for m = 1 and 
k — 3, A and 5. The innermost contour line is the interface. 




Figure 4: Validation with the exact solutions for m = 1: averaged differences 
between the numerical solution with the exact solution for /c = 3 as a function 
of the number of renormalizations n. 
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Figure 5: Exponents 5 as a function of the nonlinear exponent m, for k — 3, 
4, 5 and 6 (dots and line). The line lowermost line represents the Graveleau 
exponents P for the circular case. 
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Figure 6: Level sets of the self-similar non-circular solutions for A; = 3, 4 and 
5 for m = 1.5, m — 1.2 and m — 1.1, respectively. The innermost contour 
line is the interface. 
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Figure 7: Level sets of the initial condition (a) and an intermediate stage 
before the collapse, when the averaged radius decreased by a factor 2^ (b) 
and 2^ (c). This flow is not self-similar. 



